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Abstract. Fire emissions estimates have long been based on 
bottom-up approaches that are not only complex, but also 
fraught with compounding uncertainties. We present the de- 
velopment of a global gridded (1° x 1°) emission coeffi- 
cients (C e ) product for smoke total particulate matter (TPM) 
based on a top-down approach using coincident measure- 
ments of fire radiative power (FRP) and aerosol optical thick- 
ness (AOT) from the Moderate-resolution Imaging Spectro- 
radiometer (MODIS) sensors aboard the Terra and Aqua 
satellites. This new Fire Energetics and Emissions Research 
version 1.0 (FEER.vl) C e product has now been released 
to the community and can be obtained from http://feer.gsfc. 
nasa.gov/, along with the corresponding 1 -to- 1 mapping of 
their quality assurance (QA) flags that will enable the C e val- 
ues to be filtered by quality for use in various applications. 
The regional averages of C e values for different ecosys- 
tem types were found to be in the ranges of 16-21 gMJ -1 
for savanna and grasslands, 15-32 gMJ -1 for tropical for- 
est, 9-12 gMJ -1 for North American boreal forest, and 18- 
26 gMJ -1 for Russian boreal forest, croplands and natu- 
ral vegetation. The FEER.vl C e product was multiplied by 
time-integrated FRP data to calculate regional smoke TPM 
emissions, which were compared with equivalent emissions 
products from three existing inventories. FEER.vl showed 
higher and more reasonable smoke TPM estimates than two 
other emissions inventories that are based on bottom-up ap- 
proaches and already reported in the literature to be too low, 
but portrayed an overall reasonable agreement with another 
top-down approach. This suggests that top-down approaches 
may hold better promise and need to be further developed 
to accelerate the reduction of uncertainty associated with fire 
emissions estimation in air-quality and climate research and 


applications. Results of the analysis of FEER.vl data for 
2004-2011 show that 65-85 Tgyr -1 ofTPM is emitted glob- 
ally from open biomass burning, with a generally decreasing 
trend over this short time period. The FEER.vl C e product 
is the first global gridded product in the family of “emission 
factors”, that is based essentially on satellite measurements, 
and requires only direct satellite FRP measurements of an 
actively burning fire anywhere to evaluate its emission rate 
in near-real time, which is essential for operational activities, 
such as the monitoring and forecasting of smoke emission 
impacts on air quality. 


1 Introduction 

Smoke emitted from biomass burning is composed of a wide 
variety of particle and trace gas species that can influence 
air quality, weather, and climate variability in a significant 
way (e.g., Crutzen and Andreae, 1990; Andreae and Merlet, 
2001; Randerson et al., 2006; Schultz et ah, 2008). Among 
other sources of important atmospheric constituents (natu- 
ral and anthropogenic), open-air biomass burning is one of 
the largest contributors of both gaseous and particulate emis- 
sions to the atmosphere, and is estimated to be responsible 
for 34-38% and 40% of the global loadings of total car- 
bonaceous aerosols and black carbon (BC), respectively, as 
well as 25% of the total global carbon dioxide (CO 2 ) in- 
creases since pre-industrial times (e.g., Forster et ah, 2007). 
This is because open biomass burning occurs in most vege- 
tated parts of the world annually, in the form of natural or 
man-made burning of forests, savannas, peat lands, agricul- 
tural residues, and other ecosystem types. It is recognized 
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that an accurate understanding of smoke impacts can only be 
accomplished through accurate estimates of fire emissions 
(e.g., Langmann et ah, 2009). Therefore, researchers have 
invested considerable effort over the last several decades to 
estimate smoke emissions at different spatial and temporal 
scales from various types of biomes. Before the advent of 
satellite remote sensing, smoke emissions were estimated 
through small-scale biomass burning experiments, modeling, 
or by approximation based on proxy data such as population 
or cultural practices (e.g., Seiler and Crutzen, 1980; Hao and 
Liu, 1994; Liousse et ah, 1996). The satellite era has brought 
significant improvement in biomass burning characterization 
and emissions estimation (e.g., Dozier, 1981; Prins and Men- 
zel, 1992; Justice et ah, 1993, 2002; Cahoon et ah, 1994; 
Kaufman et ah, 1998; Giglio et ah, 2003, 2008; Wooster et 
ah, 2003, 2005; Ito and Penner, 2004; Ichoku and Kaufman, 
2005; Ichoku et ah, 2008, 2012; Schroeder et ah, 2005, 2014; 
van der Werf et ah, 2006, 2010; Giglio, 2007; Roberts and 
Wooster, 2008; Zhang et ah, 2008; Reid et ah, 2009; Vermote 
et ah, 2009; Kaiser et ah, 2012). 

Despite the considerable advancement achieved in satel- 
lite remote sensing and atmospheric modeling during the last 
couple of decades, there still remains a large uncertainty in 
the overall atmospheric impacts of aerosols and certain short- 
lived trace gases, particularly those originating from biomass 
burning such as BC and carbon monoxide (CO) (e.g., Urban- 
ski et ah, 2011; Yurganov et ah, 2011; Ichoku et ah 2012; 
Bond et ah, 2013). A major part of the uncertainty stems 
from the fact that their emission from fires are still very 
poorly constrained mainly due to the rather sporadic and 
transient character of biomass burning, which makes it diffi- 
cult to characterize experimentally (e.g., Forster et ah, 2007, 
Yokelson et ah, 2011). This can be contrasted, for instance, 
with emissions from industries and fossil-fuel burning, which 
can be quantified in a fairly straightforward manner, as the 
sources are generally stable and relatively easy to character- 
ize. For instance, the global total fossil-fuel CO 2 emissions 
are accurate to within 10% at a 95% confidence interval 
(e.g., Andres et ah, 2012), whereas the uncertainty associ- 
ated with biomass burning CO 2 emissions is still not quantifi- 
able because of lack of sufficient information (e.g., Andreae 
and Merlet, 2001), although one could infer that it would be 
as much as 100% considering the propagation of uncertain- 
ties from the various parameters that influence the emissions 
estimates (van der Werf et ah, 2010). Similar uncertainty ra- 
tios exist for other types of particulate and gaseous emissions 
from various source types (biogenic, industrial, volcanic) as 
compared to biomass burning (e.g., Diehl et ah, 2012; Bond 
et ah, 2013; Carslaw et ah, 2013). 

Many of the currently available biomass burning emis- 
sions inventories and other related products, including those 
derived from satellite data, are based on bottom-up ap- 
proaches, whereby estimates of burned biomass are derived 
from satellite-retrieved fire pixel counts, burned areas, and/or 
fire radiative power (FRP), and are then multiplied by emis- 


sion factors (EFs) of different smoke constituents derived 
from laboratory or field experiments to obtain the smoke 
emissions of these constituents (e.g., Chin et ah, 2002, Ito 
and Penner, 2004; Hoelzemann et ah, 2004, Liousse et ah, 
2004; Michel et ah, 2005; van der Werf et ah, 2006, 2010; 
Generoso et ah, 2007). Examples of such inventories that are 
currently being used by the community include GFED (van 
der Werf et ah, 2006, 2010), GFAS (Kaiser et ah, 2012), 
FLAMBE (Reid et ah, 2009), FINN (Weidinmyer et ah, 
2011), and GBBEP-Geo (Zhang et ah, 2012). Recent re- 
search findings suggest that such bottom-up approaches lead 
to severe underestimations particularly of smoke aerosols un- 
less bias correction is applied through modeling (e.g., Li- 
ousse et ah, 2010; Kaiser et ah, 2012). Top-down approaches 
are starting to be investigated for deriving biomass-burning 
emissions, sometimes in conjunction with model assimila- 
tion (e.g., Sofiev et ah, 2009; Kaiser et ah, 2012; Darmenov 
and da Silva, 2013). Although biomass burning emits several 
dozens of particulate and gaseous species (e.g., Andreae and 
Merlet, 2001; Akagi et ah, 2011), this study is specifically 
focused on smoke aerosol or total particulate matter (TPM) 
emissions. 

This paper presents the development of the first gridded 
global top-down biomass burning aerosol emission coeffi- 
cients product that is based strictly on locally collocated 
satellite measurements of both fire radiative power (FRP) 
and aerosol optical thickness (AOT). The original idea and 
an initial algorithm were developed in Ichoku and Kauf- 
man (2005) in which FRP and AOT retrieved from the 
Moderate-resolution Imaging Spectro-radiometer (MODIS) 
sensor aboard the NASA Terra and Aqua satellites were uti- 
lized together with wind vectors from the National Cen- 
ter for Environmental Prediction/National Center for Atmo- 
spheric Research (NCEP/NCAR) meteorological reanalysis 
data to generate smoke-aerosol emission coefficients (C e in 
kgMJ -1 ) for several biomass burning regions. Such top- 
down emission coefficients are found to be useful, as simply 
multiplying C e by the satellite-retrieved FRP of a fire gives 
the corresponding instantaneous TPM emission rate for that 
fire. Likewise, in the case of consistent and frequent fire ob- 
servations such as from a geostationary platform, multiply- 
ing C e by the time-integrated FRP (or fire radiative energy, 
FRE) gives the TPM emission for that time interval. This 
C e x FRP (or C e x FRE) emissions estimation approach and 
variants of it have been subsequently developed and imple- 
mented successfully in various regional studies (e.g., Jor- 
dan et ah, 2008; Henderson et ah, 2008, 2010; Sofiev et al„ 
2009; Vermote et ah, 2009). However, the original Ichoku 
and Kaufman (2005) algorithm has been substantially en- 
hanced and used to generate a gridded global C e product 
using an updated algorithm and newer versions of FRP and 
AOT data from MODIS as well as wind vectors from the 
Modem Era Retrospective- Analysis for Research and Appli- 
cations (MERRA) data sets (Rienecker et ah, 2011) provided 
by the NASA Goddard Global Modeling and Assimilation 
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Office (GMAO). The newly generated gridded C e data prod- 
ucts are available at the NASA Fire Energetics and Emis- 
sions Research (FEER) web site (http://feer.gsfc.nasa.gov/) 
together with MODIS FRP data and links to other relevant 
satellite FRP data. 

Section 2 provides the background and theoretical basis 
of the approach. Section 3 describes the characteristics of 
the various input data (FRP, AOT, winds) used to calculate 
C e . Section 4 gives the full details of the updated method- 
ology for deriving C e and the associated uncertainty analy- 
ses. Section 5 presents the use of the gridded C e product to 
estimate smoke particulate emissions over different regions 
and comparisons with similar emission products, namely, 
the Global Fire Emissions Database version 3.1 (GFED.v3: 
van der Werf et al., 2006, 2010), the Global Fire Assimila- 
tion System version 1.0 (GFAS.vl: Kaiser et ah, 2012), and 
the Quick Fire Emission Dataset version 2.4 (QFED.v2: van 
Donkelaar et ah, 2011; Darmenov and da Silva, 2013). Fi- 
nally, Sect. 6 provides a summary and relevant concluding 
statements. 

2 Background and theoretical considerations 

Traditionally, the amount of a given aerosol or trace-gas 
species emitted from open biomass burning is derived by 
multiplying that species’ emission factor (in grams of species 
per kilogram of dry matter burned) by the mass of biomass 
burned. The basic equation is of the form (e.g., Andreae and 
Merlet, 2001): 

M x — EF X • Mbiomass . (1) 

where M x is the mass of emitted smoke species x, EF X is the 
emission factor for the emitted species x, and Mbiomass is the 
mass of the dry biomass burned. 

A similar relationship to Eq. (1) was established by Ichoku 
and Kaufman (2005) in which EF X is replaced with C* , 
which is designated as the emission coefficient (for any given 
species x), and Mbiomass is replaced with either FRE or its re- 
lease rate Rf re (i.e., FRP). Thus, 

M x = C* ■ FRE 

or 

//x = c e v • R fre , (2) 

where R x is the rate of emission of species x (expressed 
in kg s 1 ) since Rf re is the FRE release rate expressed in 
MJs -1 , or MW. Cg is therefore expressed in kgMJ -1 . The 
validity of the relationship in Eq. (2) has been verified in 
a laboratory experiment, where satellite measurements of 
fire energetics and smoke were replicated by burning small 
biomass fuel samples in a burn chamber equipped with a gi- 
ant smoke stack upon which the relevant instruments were 


set up, and the retrieved FRP and AOT were used to derive 
C e for smoke aerosols (Ichoku et al., 2008b). 

Equations (1) and (2) are functionally very similar, and re- 
lating the two would suggest that there is a linear relationship 
between Mbiomass and FRE. Indeed, a series of field experi- 
ments showed that FRE is proportional to Mbiomass in a linear 
fashion, such that Mbiomass = 0.368(± 0.015)- FRE, in which 
the numeric coefficient (0.368 kg MJ -1 ) is designated as the 
biomass consumption factor (F c ) (Wooster et al., 2005). The 
Wooster et al. (2005) study indicated that the same relation- 
ship is expected to hold for satellite observations when total 
biomass consumed Mbiomass is substituted with the rate of 
biomass consumption and FRE with Rf vc . That relationship 
has also been verified in laboratory experiments (Freeborn 
et al., 2008; Ichoku et al., 2008b), and has been applied in 
the estimation of Mbiomass over Africa using FRE derived by 
integrating Rf re measurements from the Spinning Enhanced 
Visible and Infrared Imager (SEVIR1) aboard the Meteosat 
second generation (MSG) series of European geostationary 
meteorological satellites (Roberts et al. 2005, 2011). Sim- 
ilarly, as derived in Ichoku et al. (2008b), the mass-based 
emission factor, EF X , in Eq. (1) is related to the FRE-based 
emission coefficient C) for any given fire-emitted species x 
as 

EF X = Cg / F c , (3) 

where F c is the biomass consumption factor ( F c ) defined in 
Wooster et al. (2005). 

This ability to relate the satellite-measured fire radiant 
heat release rate Rf re and the top-down derived emission 
coefficient C e to physical quantities of combusted biomass 
Mbiomass and its associated bottom-up smoke emission factor 
EF X , respectively, is a major motivation buttressing the study 
described in this paper. Currently, only a few generalized val- 
ues of EF X are available for certain ecosystem types, which 
is highly limiting given that EF X is likely to vary by location 
in the same manner as fuel characteristics, even within the 
same ecosystem type. Therefore, by using satellite-measured 
Rf re and smoke aerosols to derive C e globally as a gridded 
product based on the developed top-down approach, it is not 
only possible to compare these results with those based on 
bottom-up approaches, but it can even lead to the develop- 
ment of a gridded EF X product that would offer a much finer 
spatial coverage and resolution than do the current products. 

3 Data 

The main data products used in generating the gridded C e 
are satellite measurements of FRP and AOT, as well as as- 
similated wind fields from MERRA. Both the FRP and AOT 
products used in this work are derived from the MODIS sen- 
sors aboard the (1) Terra satellite launched in 1999 with local 
equator crossing times of 10:30 a.m. and 10:30 p.m., and (2) 
Aqua satellite launched in 2002 with local equator crossing 
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times of 1:30 p.m. and 1:30 a.m. The analysis in this paper is 
based on data covering the period between 2003 and 2010, 
inclusive. The specific attributes of these products, such as 
their spatial and temporal resolutions, versions, and uncer- 
tainties are discussed in the following sub-sections. It should 
be noted that MOD1S data versions are essentially referred 
to as data “collections”, a terminology that will be used 
throughout this paper. 

3.1 Fire radiative power 

Active fire observation products from MOD1S on Terra 
(MOD14) and Aqua (MYD14) are provided at a nominal 
spatial resolution of 1 km at nadir (Justice et al., 2002; Giglio 
et ah, 2003). FRP (or Rf re ) is one of the main parameters pro- 
vided within these products for every fire pixel detected. The 
original formulation for derivation of Rf le was developed in 
Kaufman et ah (1998, p. 32226, Eq. (1)) as 

tffre = 4.34 x 10~ 19 • (r 4 8 - r 4 8 b , ) (4) 

where R& e is the pixel fire radiative power (in Wm -2 ), T 4 
is the fire pixel brightness temperature (in K) at the 4 pm 
channel (3.96 pm for MOD1S), and T 41 , is the 4 pm brightness 
temperature of the background surrounding the fire pixel. 

Equation (3) was used to derive FRP values from MOD1S 
up to the collection 4 data set released in 2004. Those collec- 
tion 4 data were used for the Ichoku and Kaufman (2005) 
study. Starting from collection 5, the right hand side of 
Eq. (3) was multiplied by the area of each pixel to account for 
the variation of ground pixel size with MODIS scan angle, re- 
sulting in units ofMW for Rf re (Giglio, 2010). The collection 
5 FRP data set, which is the latest data version available at the 
time of this study, has been used for the calculations reported 
here. The potential effects that this change in FRP values has 
on computed C e is analyzed in Sect. 4.6. However, it is note- 
worthy that FRP retrievals from MODIS have not yet been 
validated, even though the uncertainty associated with the 
detection of fire locations has been characterized using fire 
detections at 30 m nominal spatial resolution from the En- 
hanced Thematic Mapper Plus (ETM+) sensor aboard the 
Landsat-7 satellite and the Advanced Spaceborne Thermal 
Emission and Reflection Radiometer (ASTER) aboard Terra 
(e.g., Morisette et al., 2005a, b; Schroeder et al., 2008a, b). 

3.2 Aerosol optical thickness 

The AOT (r a x) data used for this study were also re- 
trieved from MODIS on Tern (MOD04 L2) and Aqua 
(MYD04 L2) at 10 km spatial resolution at nadir. MODIS 
measures AOT at 470, 550, 660, and 2100nm wavelengths 
(X) over land, and at 470, 550, 660, 870, 1200, 1600, and 
2100 nm wavelengths over ocean (e.g., Remer et al., 2005, 
2008; Ichoku et al., 2005; Levy et al., 2010). However, only 
the AOT data retrieved over land are used in this study, 
since smoke from biomass burning can only be emitted over 


land where fires normally occur, although this makes it dif- 
ficult to get a sufficient amount of retrievals for fires oc- 
curring very close to coastlines. Specifically, we use AOT 
measurements at 550 nm wavelength, as this falls within the 
mid-visible or green region of the electromagnetic spectrum, 
which is the most commonly used wavelength region in 
aerosol radiation studies. Unlike the FRP data, MODIS AOT 
data have been extensively characterized and validated using 
ground-based sun-photometer measurements from the global 
Aerosol Robotic Network (AERONET, e.g., Holben et al., 
1998, 2001). However, as with the fire product, the collection 
5 MODIS Level 2 Aerosol Product (http://modis-atmos.gsfc. 
nasa.gov/C005_Changes/C005_Aerosol_5. 2.pdf) was used 
in this study instead of the collection 4 that was used in 
Ichoku and Kaufman (2005). Detailed analyses of the effects 
of this change in AOT data version and other factors on the 
computed C e results are presented in Sect. 4.6. 

3.3 Wind vectors 

The wind vectors used for this study were extracted from 
MERRA’s inst3_3d_asm_Cp product provided at a spatial 
resolution of 1.25° x 1.25° and a temporal resolution of 
3 h. The documentation for that product is available at 
http://disc.sci.gsfc.nasa.gov/mdisc/data-holdings/merra/ 
inst3_3d_asm_Cp.shtml. Wind data at pressure levels of 
925, 850, and 700 mb, roughly corresponding to heights 
above mean sea level (a.s.l.) of 750 m, 1.5 km, and 3 km, 
respectively, were extracted and used for spatial aerosol data 
analysis to derive smoke TPM emission rates. However, after 
the analyses, the wind data at 850 mb (roughly 1.5 km a.s.l.) 
were used to generate the final product as described in 
Sect. 4 and in Ichoku and Kaufman (2005), because most 
fires observable from satellite at ~1 km spatial resolution 
in different parts of the world typically inject plumes to 
heights of about 1.5±1. 0km above ground level (a.g.l.), with 
the exception of very large fires (e.g., Lavoue et al., 2000; 
Freitas et al., 2006; Kahn et al., 2007; Val Martin et al., 
2012; Yang et al., 2013). 

3.4 Other data 

Several other data types, products, and parameters were used 
in this study. The global average aerosol mass extinction ef- 
ficiency value of j 8 e = 4.6 m 2 g _ 1 that was used in Ichoku and 
Kaufman (2005), based on the work of Reid et al. (2005), has 
also been used in the current work. Coincident digital eleva- 
tion model (DEM) data and ecosystem data were obtained 
for each data point during processing for reference. DEM 
data at 30arcsec resolution (GTOPO30: https://lta.cr.usgs. 
gov/GTOPO30) are provided by the US Geological Surveys 
(USGS). We used the DEM data sets to determine the land- 
surface elevation, land-sea mask, slope, and aspect. Ecosys- 
tem data used in this work are from the 1 arcmin resolution 
global ecosystem map of 2004 derived from MODIS (http: 
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//modis-atmos. gsfc.nasa.gov/ECOSYSTEM/) using the In- 
ternational Geosphere/Biosphere Program (IGBP) classifica- 
tion scheme. Digitized smoke plume data from the Multi- 
angle Imaging Spetro-Radiometer (MISR) INteractive eX- 
plorer (MINX) tool (Nelson et al., 2008, 2013) were used 
to evaluate the relationship between the wind direction from 
MERRA and the actual plume direction as observed on the 
MOD1S imagery. 


4 Methodology 

The basic methodology for deriving the smoke-aerosol emis- 
sion coefficients C e from satellite measurements of Rf re and 
t a x was developed in Ichoku and Kaufman (2005). However, 
although the basic structure and processing sequence of the 
original algorithm have been maintained, several adjustments 
and updates were required, in terms of both the algorithm and 
input data, in order to generate the gridded products reported 
in this paper. 

4.1 Algorithm logic for C e calculation 

The logic progression within the algorithm to calculate C e is 
generally similar to that described in Ichoku and Kaufman 
(2005) in that the first stage of the algorithm is completed 
on a 10 km resolution aerosol-pixel level, followed by a sec- 
ond stage where these units are aggregated within larger ar- 
eas (1° x 1° regular grid in the present case), and then ending 
with the actual calculation of C e . The core of the algorithm is 
outlined in this section and visualized in Fig. 1. The specific 
details of the three analysis stages, including the main adjust- 
ments and updates applied in the current study, are described 
in Sects. 4.2, 4.3, and 4.4. 

The first stage of the algorithm is designed to generate val- 
ues of Rf re and R sa (the smoke-aerosol emission rate) for 
each aerosol pixel with detected fire(s). Fitting the MODIS 
1 km resolution active fire data into the corresponding 10 km 
resolution aerosol-pixel data is very straightforward because 
both data sets originate from the same instrument on the same 
platform and from the same original data product. Therefore, 
the Rf rc for a given aerosol pixel is derived as 


Nf 

R fr e = X>RP ; , (5) 

(=1 

where FRP is the fire radiative power measurement of indi- 
vidual active fire pixels, and Nf is the total number of active 
fire detections within a given aerosol pixel. 

Derivation of R sa is less straightforward and involves cal- 
culations utilizing AOT and wind vectors in a 3 x 3 aerosol 
matrix centered on the fire-affected aerosol pixel, as depicted 
in Fig. 2. Since the plume can easily influence neighbor- 
ing pixels, the 3x3 matrix is split into four quadrants, one 
of which is deemed to be the “downwind quadrant” based 
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Figure 1. Flowchart of the core algorithm for deriving FEER.vl 
emission coefficients (C e ) of smoke total particulate matter (TPM), 
as outlined in Sect. 4.1, showing only the relevant input and output 
variables/parameters and their associated equation numbers. This 
does not include the various sensitivity studies, data selection, re- 
finement, and filtration processes, as well as the post-processing 
gap-filling and evaluation steps that are described in the rest of 
Sect. 4. 


on the wind direction, and the four pixels therein are as- 
sumed to contain parts of the plume, whereas the remain- 
ing five pixels are the background. The zonal (u) and merid- 
ional (i>) components of the wind vector at the 850 mb atmo- 
spheric pressure level are used to calculate the wind speed 
(WS = Vm 2 + v 2 ) and azimuth (0 = cos -1 [u/WS]). The 
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Figure 2. Spatial configuration of a 3 x 3 aerosol-pixel matrix lay- 
out, whose central pixel contains fires, showing the four downwind 
pixels (shaded red, quadrant IV) classified as having smoke, and the 
five remaining upwind pixels (shaded blue) constituting the back- 
ground. The downwind quadrant is determined by the wind direc- 
tion. The pixel indices (0-8) shown in their bottom left-hand cor- 
ners are defined by their scanning configuration, signified here by 
the directions of line and sample coordinates. The sample direction 
is along-scan and the line direction is along-track. (The background 
image taken by Aqua/MODIS at 20:45 UTC on 1 July 2012 shows 
the Fontenelle Fire in Wyoming, USA.) 


azimuth is compared with the MODIS along-track direction 
to determine the downwind quadrant in which the plume is 
located (see Fig. 2). The AOT that is attributable to the fire(s) 
within the central aerosol pixel can be determined as 


l a550 


= r. 


a550 ' 


l a550’ 


( 6 ) 


where the superscripts f, t, and b, respectively, designate the 
fire-emitted, total, and background AOT at 550 nm wave- 
length. The background AOT value, r a b 550 , is calculated as 
the mean of the valid background AOT values (shown in blue 
in Fig. 2), weighted by aerosol-pixel area. The fire-emitted 
AOT, r b 550 , is found by subtracting this mean r b 550 value 
from r* 550 of each aerosol pixel in the downwind quadrant 
(plume direction). Individual r a550 values that are negative 
are set to zero. Subsequently, an area-weighted average r,| 55(| 
value is calculated from the downwind aerosol pixels to rep- 
resent the unit plume being analyzed. Thus, 



where A is the aerosol-pixel area, N a t is the number of valid 
aerosol retrievals in the downwind quadrant, and N a f is the 
number of valid aerosol retrievals in the downwind quadrant 


whose r* 550 appropriately exceeds r b 550 . This fire-emitted 
AOT ( r.[- 5() ) is converted to smoke -aerosol column mass den- 
sity (Md in gm -2 ) as 

M d = T a f 550/ / g e , ( 8 ) 

where /3 e (expressed in m 2 g -1 ) is the smoke aerosol specific 
extinction or mass extinction efficiency derived from Reid et 
al. (2005). Using the total area of the four downwind pixels. 
At , the mass of smoke-aerosol emission is then calculated 
by 

M sa — Md ■ A t . (9) 

Determining the smoke-aerosol emission rate R sa requires 
knowledge of how much time, T, it must have taken to emit 
Md- For a given plume, T is assumed to be the time it would 
take for the wind to clear all smoke aerosol from the down- 
wind quadrant within the 3x3 aerosol-pixel matrix, and is 
estimated as 

T = L /ws> ( 10 ) 

where L represents the length of the plume within the 3x3 
aerosol-pixel matrix. In the case where there are multiple ac- 
tive fire detections within one aerosol pixel, the plume dis- 
tances are averaged to yield one value for L. Finally, the rate 
of smoke-aerosol emission is estimated as 

Rs a = Msa /:T’ (H) 

where R sa is expressed in kg s _1 . Thus, paired values of R r re 
and R sa are calculated for each fire-containing aerosol pixel 
and collected into a “pixel-level” product, which is used in 
the second stage to generate similar measurements at larger 
scales. 

The second stage of the algorithm involves aggregating the 
pixel-level calculations of Rf re and R sa to determine corre- 
sponding values for larger areas or regions at each MODIS 
overpass event. Since it is the aim of this study to render the 
FEER.vl C e product in a 1° x 1° - grid configuration, corre- 
sponding pixel-level Rf re and R sa values within each 1° x 1° 
grid cell are aggregated as 

*fie = E/4e ( 12 ) 

and, 

r™=J 2, r ™- ( 13 ) 

Details of the implementation of this stage, including its 
unique sensitivity analyses and data filtering processes, are 
described in Sect. 4.3. The output is the “grid-level” data set 
to be used for the third and final stage during which C e is 
actually derived. 

The procedure to derive C e is to generate a scatter plot 
for each 1° x 1° grid cell using the Rf re and R sa data calcu- 
lated in the second stage for a specified time domain (in our 
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(a) FRP [MW] (b) FRP (MW] 

Figure 3. Scatter plots of smoke emission rate (R sa ) against fire radiative power (FRP or Rf re ) derived from both Terra and Aqua MODIS 
observations during the period 2003-2010 for a 1° x 1° grid cell centered at (a) 1.5° S, 15.5° E, and (b) 22.5° N, 115.5°E. The latter 
demonstrates the effect of removing outliers in such scatter plots. The outlier is identified in red and the blue line is the linear least-squares 
regression fit through the remaining points, which in this case results in a C e value of 0.0747 and an r~ value of 0.82. This is a great increase 
over the case without the outlier removal process, whose regression line is shown in gray and has much lower values of both C e (0.0128) and 
level of confidence (r 2 = 0.16). 


case 2003-2010), with R sa on the dependent axis. A mini- 
mum of six points is allowed for a scatter plot, as we con- 
sider six to be a minimally reasonable number of data points 
for linear regression fitting. A zero-intercept regression line 
(of the form y =mx, where y — R sa and x — R[ IC ) is fitted to 
the scatter plot for each grid cell (Fig. 3). The gradient m is 
the coefficient of emission C e and the goodness of fit is eval- 
uated on the basis of the coefficient of determination (r 2 ), 
as will be described in Sect. 4.4. Hence, for grid cells with 
good fits, Rf re only needs to be multiplied by C e to derive the 
smoke emission rate R sa , even in near-real time, bearing in 
mind the possibility of large biases due to the inherent differ- 
ences in individual fire characteristics even within the same 
fire regime (e.g., Schroeder et al., 2014). 

4.2 First stage: pixel-level data analysis 

The calculations of Rf re and R sa begin at the pixel-level. 
Several changes have been made in the current study rela- 
tive to the original method described in Ichoku and Kaufman 
(2005). In that paper, r a550 was defined as the maximum AOT 
measured in the 3x3 matrix, whereas r b 550 was defined as 
the minimum AOT measured in the eight pixels immediately 
surrounding the center pixel, regardless of the actual direc- 
tion of the plume. That methodology should produce good 
results when the plume is prominent and the background is 
uniform and clear. Otherwise, such as when the plume is 
thin or highly dispersed, or when plumes from a different 
fire enter any of the aerosol pixels within the 3x3 matrix, 
the result can be unpredictable. To characterize such situa- 
tions, 240 digitized MISR plumes from fires that occurred 
in Siberia in May 2003 were analyzed. The outlines of the 
3x3 matrix of MODIS aerosol pixels centered on each fire 
were delineated over a MISR true-color imagery, and the cor- 


responding MODIS AOT values were recorded along with 
wind directions from both MISR (observed) and MERRA 
(modeled). For each of the 240 surveyed plumes, visual clas- 
sifications were made with the help of the MISR fine (275 m) 
spatial-resolution imagery to identify which 3x3 matrices 
of MODIS aerosol pixels contained: plumes, clouds, haze, or 
fires. Of the cases analyzed, 64 % had detected fires in the 
surrounding pixels and 61 % had background smoke or haze. 
These proportions can be quite different in other regions be- 
cause fire density and smoke dispersion characteristics vary 
by region, biome, and season. In particular, regions that typi- 
cally have relatively smaller fire sizes, such as the African sa- 
vannas, are likely to be more impacted by this phenomenon. 
This significant percentage of extraneous aerosol contami- 
nation can have an adverse impact on the determination of 
r a550 alu * r a550’ anc * consequently also on the accuracy of 
r a550 ( see 7). For instance, in the Ichoku and Kaufman 
(2005) method: (i) if a larger plume exists in a neighboring 
pixel its AOT can be erroneously taken as the r* 550 for the 
current plume, and (ii) unless any external smoke entering 
into the 3x3 matrix affects every one of the aerosol pixels 
similarly, the minimum AOT eventually used as r b 550 may be 
much lower than the average background AOT. Furthermore, 
in case (ii) r a550 may be likely boosted to a higher value, po- 
tentially resulting in the overestimation of r a550 when r b 550 
is subtracted from r.j 55(| . To mitigate such situations in the 
current study, the plume direction is first identified based on 
the relative wind direction, as shown in Fig. 2, whereupon a 
combination of the four downwind pixels are used to calcu- 
late more realistic values of r* 550 . Likewise, instead of using 
the overall minimum AOT, the average AOT of the five up- 
wind pixels are used to more realistically determine r a b 55() . It 
is assumed that where smaller smoke plumes from fires in the 
neighboring or external pixels are present, the distribution of 
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their effects among the four downwind and/or five upwind 
aerosol pixels will tend to dilute the smoke contamination on 
T a550 an d r a550 va l ues > thereby minimizing its impact on the 
resulting r,| 5 - (| value. Cases that have more serious contami- 
nation also have a good chance of being eliminated both by 
the threshold requirements outlined in Sect. 4.3 and the out- 
lier removal process described in Sect. 4.4. However, another 
issue that can potentially affect the result, though it is not di- 
rectly related to this t,[ 55() algorithm, is the fact that near- 
source thick smoke plumes are occasionally misclassified 
as clouds and therefore omitted in the collection 5 MOD1S 
aerosol product (e.g., Livingston et al., 2014; Schroeder et 
ah, 2014). For instance, if any of the four downwind pixels 
in a 3 x 3 matrix has no AOT value, this could lead to un- 
derestimation of r* 550 and consequently also of M sa and R sa 
for the specific plume unit being analyzed. Such situations 
are expected to improve with future enhancements in aerosol 
retrieval algorithms particularly close to smoke sources, but 
the filtering and outlier removal processes implemented in 
this study are helpful in the meantime. 

One important condition in classifying downwind and up- 
wind sections is that the wind direction needs to be correct. 
The level of accuracy, however, is variable since the actual 
requirement is that only the correct downwind quadrant is 
identified. The MINX data set for Siberia in May 2003 also 
makes the evaluation of this condition possible, and showed 
that the success rate of using MERRA to correctly identify 
the downwind quadrant was on the order of 80 %. This is an 
acceptable rate, especially considering the fact that most of 
the failed cases will likely be filtered out in the second stage 
of the C e algorithm (Sect. 4.3) due to a probable decrease in 

T a550 anc * i ncrease i' 1 r a550 suc ^ ^ at T a550 W 'H to ° l° w - ^ 
should also be noted that there was no increase in accuracy 
when data were matched to the closest plume injection level 
as recorded in the MINX database than when only the 850 mb 
pressure level data were used. This reaffirms the validity of 
the use of wind data at 850 mb for generating the FEER.vl 
C e product, although since this is based only on data from 
Siberia, it may not be ideal for other parts of the world where 
smoke plumes are typically injected either much lower or 
much higher than the 850 mb pressure level. Examples in- 
clude the African Sahel region where fires are typically small 
in sizes and inject plumes lower than 1 km (e.g., Yang et al., 
2013) and the Canadian boreal forest region where very large 
fires can inject plumes higher than 2 km (e.g., Lavoue et al., 
2000; Val Martin et al., 2012). 

Another measure taken to minimize uncertainty in the 
pixel-level analysis relative to the original algorithm in 
Ichoku and Kaufman (2005) is the use of the wind vectors in 
the derivation of the distance ( L ) the plume travels within the 
analysis unit (the 3x3 aerosol-pixel matrix). The original 
method equates L to the square root of the central aerosol- 
pixel area without considering the actual relative positions 
of the individual 1 km resolution fire pixels within the 10 km 



0 30 60 90 120 150 180 210 

Plume Time [minutes] 

Figure 4. Probability density functions of plume time (Eq. 10) for 
different filter configurations (Table 2) used in data screening and 
selection for generating C e scatter plots. 

resolution aerosol pixel. The new algorithm takes into ac- 
count the relative positions of these fire pixels within the cen- 
tral aerosol pixel in estimating the distance traveled by each 
smoke plume from its source (center of the fire pixel) to the 
edge of the 3x3 aerosol-pixel matrix (see Fig. 2). L is ex- 
tended to the edge of the 3x3 pixel matrix, instead of only 
to the edge of the central aerosol pixel, to prevent any ambi- 
guity in L from introducing large errors in R sa calculations 
(Eqs. 10 and 11), particularly when the fire is very close to 
the downwind edge of the central aerosol pixel. Therefore, 
provided the smoke plume follows MERRA’s wind direc- 
tion at 850 mb, it is believed that the derived values for L 
and consequently T, R sa , and C e will be much more accu- 
rate. Although conceptually this algorithm change is an im- 
portant improvement, however, for relatively small fires and 
low wind-speed situations, plumes may not reach the edge of 
the 3x3 aerosol-pixel matrix, resulting in overestimation of 
L and consequently T (Eq. 10). Figure 4 shows the proba- 
bility density functions (PDF) of the plume time, T, for the 
four data-filter settings (00000, 10000, 11000, 11300: low- 
est to highest quality) used in this study, as explained later 
in Sect. 4.3. Overall, the maximum probability for T lies in 
the range of 45-55 min, with a gradual decrease beyond the 
peak. For such large time periods of R sa estimation, the fire 
that emitted r,|- 5(| may have changed significantly relative to 
the FRP value recorded at the time of observation. This is the 
weakness of using the instantaneous and simultaneous ob- 
servations of a fire and corresponding plume based on a rela- 
tively low spatial-resolution aerosol product. However, if the 
spatial resolution of the input aerosol product improves (as is 
currently being developed by the MODIS aerosol team), this 
issue will be alleviated. 

Lastly, in the original algorithm by Ichoku and Kaufman 
(2005), single values of R sa and Rf Te were calculated for large 
regions or areas (in this case 1° x 1° grid cells) involving 
multiple fire/plume units only after the upstream variables 
had been aggregated into these regions. That approach has 
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Figure 5. Selected 1° x 1° grid cells fora sensitivity analysis on C e 
scatter plots and values based on using different threshold param- 
eters and settings are identified on this MODIS true-color image. 
These sites were selected with the intention of maintaining diversity 
in location, fire type, biome, number of data points, and expected 
goodness-of-fit of linear regression. 


been modified in the current implementation to minimize its 
vulnerability to errors that may be inherent in the aggrega- 
tion processes preceding the calculations. In the current al- 
gorithm, the pixel-level analysis is continued up until the cal- 
culation of R sa and A’f rc for each fire/plume unit. This allows 
for flexibility in the use and aggregation of these products at 
different scales and corresponding uncertainty estimation. In 
the current work, the values of R sa and Rf le generated at the 
pixel-level are aggregated into the 1° resolution grid cells for 
creating scatter plots. 

4.3 Second stage: gridded data analysis 

The creation of a gridded product at 1° resolution arises from 
the need for derivation of a gridded smoke emission coeffi- 
cient C e that would be available for use in generating emis- 
sions wherever fires occur around the world for various types 
of applications and modeling. Because the pixel-level smoke- 
aerosol emission rates parameter (R sa ) simply reports values 
for all aerosol pixels containing fire regardless of the qual- 
ity of the aerosol retrievals, the development of this grid- 
ded product necessitates a methodology for removing invalid 
or erroneous data, which is accomplished through the use 
of thresholds applied to selected parameters. These are de- 
scribed in Table 1, along with the purpose for using each one 
of them. 

To determine appropriate thresholds for these parameters, 
several 1° x 1° grid cells were selected around the globe 
from a variety of biomass burning regions to conduct sensi- 
tivity analyses using data from the full time period of 2003- 
2010 (Fig. 5). The selections were made by examining ran- 
dom grid cells spread out throughout the entire globe and 
manually ensuring that the final selection maintained diver- 
sity in location, fire type, biome, number of data points, 


and expected goodness-of-fit of linear regression. Data con- 
tained within these sample grid cells were used to perform 
a dynamic, detailed analysis of the calculations described in 
Sect. 4.1 (and illustrated in Fig. 1) to quickly generate differ- 
ent emission coefficients. For each site, these algorithmic cal- 
culations to aggregate pixel-level values of R& e and R sa into 
the grid cell and to calculate C e were applied inside an Excel 
workbook, where provisions were made for a user to control 
the threshold parameters listed in Table 1 . Each threshold pa- 
rameter was varied and studied in different combinations as 
their effects on the final results were visualized. The calcula- 
tions were followed through all the way to the scatter plots of 
R sa and Rf ie , and a linear least-squares regression line pass- 
ing through the origin was fitted, resulting in values of C e . In 
this way, the corresponding change in the look of the scatter 
plot and in the value for C e due to varying threshold settings 
was observed in real time. Thus, the results were dynamic in 
nature and allowed for proficient sensitivity analysis at each 
of the sites. 

A five-digit code was developed to represent the different 
combinations of the threshold settings, as designated in the 
header row of Table 2. Each digit within the five-digit code 
represents one set of parameters that are changed, and the 
digit number represents different settings for those param- 
eters. Thus, the 00000 setting represents the case when no 
filtering is applied to the data set at all, except the standard 
requirement that there be valid retrievals of Rf re (F_power) 
and R sa . A basic set of parameters were selected as a com- 
mon improvement in all the selected sites, identified as a “1” 
for the first digit in the settings code (i.e., 10000 is a set- 
ting with only these basic settings turned on). These basic 
parameters are (Tables 1 and 2) the scan angle of the aerosol 
pixel, the wind speed, the number of available surrounding 
aerosol retrievals, the lowest AOT quality assurance flags 
specified for selecting the upwind and downwind aerosol pix- 
els, and the number of valid AOT retrievals of the upwind 
and downwind pixels. The second digit of value “1” (i.e., 
11000) represents the elimination of cloud contamination by 
setting A cloud fraction mean to 0. This setting produced 
the largest single noticeable improvement across the board, 
not only in reduced point scatter, but also in improved re- 
gression line fits. The third digit setting corresponds to the 
next set of thresholds used to impose restrictions on extreme 
minima in the main parameters contributing to the calcula- 
tion of C e , namely r a550 (A_AOT550_fire) and Rf Ie . Over 
the course of examining sufficient threshold values to use 
for these parameters, two values for each parameter were 
selected for further testing with all the sites collectively, 
creating four possible combinations: “1” (r^ 550 >0.01 and 
Rfre > 15 MW), “2” (r a f 550 > 0.01 and fl fre > 20 MW), “3” 
(r a f 550 > 0.02 and R fle > 15 MW), and “4” (r a f 550 > 0.02 and 
Rfre > 20 MW). This was motivated by the realization that 
extremely low r a550 and Rf re values within a 10 km x 10 km 
aerosol pixel would be too close to the noise level to be good 
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Table 1 . List of parameters that are used for data filtering in the gridded product development step described in Sect. 4.3 (parameter-name 
prefixes “A”, “F” and “M” indicate whether a parameter belongs to the MODIS Aerosol, MODIS Fire, or MERRA Meteorological data sets, 
respectively.) 


Parameter 

Description 

Purpose 

A scan angle 

The scan angle of the aerosol pixel. 

Eliminate the effect of pixel overlap, which adds too 
much complexity in determining total and upwind AOT 
values. 

M wind speed 

The wind speed from MERRA. 

Eliminate slow air mass that would escalate T and make 
Rsa very small. 

A retrievals nearby 

The number of available aerosol retrievals immediately 
surrounding the center pixel. 

Ensure the pixel is not along the edge of the MODIS 
scene (granule), and that no nearby feature prevents 
aerosol retrieval. 

A_AOT550_fire 

Fire-emitted AOT at 550 nm (i.e., total - background 
AOT at 550 nm). 

Eliminate cases where the plume signal is weak relative 
to the background. 

A QA AOT total 

The smallest QA used in selecting total AOT from the 
downwind pixels. 

Allow flexibility to specify desired range of AOT quality 
flags. 

A_QA_AOT bkgd 

The smallest QA used in selecting background AOT 
from the upwind pixels. 

Reduce uncertainty in background (upwind) AOT mea- 
surements. 

A_AOT550_retr_total 

The number of valid downwind aerosol retrievals. 

Ensure that enough valid pixels are available for accu- 
rate total AOT determination. 

A_AOT5 50_retr_bkgd 

The number of valid upwind aerosol retrievals. 

Ensure that enough valid pixels are available for accu- 
rate background AOT determination. 

A cloud frac mean 

Mean cloud fraction of the 3x3 aerosol-pixel matrix 
for unit plume analysis. 

Reduce the chances of cloud being falsely identified as 
smoke and vice versa, or cloud obscuration of fire. 

F_pcounts 

The number of MODIS fire pixels inside the center 
aerosol pixel. 

Optimize number of fire pixels within aerosol pixel for 
accurate FRP total. 

Fjcounts nearby 

The number of MODIS fire pixels in surrounding 8 
aerosol pixels. 

Eliminate uneven contamination of AOT by emissions 
from nearby fires. 

F_pcounts DW3 

The number of fire pixels in three downwind pixels (ex- 
cluding center). 

Eliminate contamination of target plume by those from 
nearby fires. 

Fjiower 

The cumulative FRP value of all fires within the center 
aerosol pixel. 

Limit small fires and underestimated FRP values that 
can cause large errors. 

Rsa 

The rate of smoke emission. 

Limit invalid values or cases with insignificant amounts 
of smoke production 


for useful analysis. However, between the two values of the 
T a 550 threshold tested, 0.02 was adopted as more realistic for 
further analysis because it is closer to the absolute compo- 
nent (i.e., 0.05) of the expected AOT retrieval error ± (0.05 
+ 15%) from MODIS over land (e.g., Levy et al., 2010). 
Also, by observing the effect of different choices of /Are 
thresholds on the sites collectively, it became visually appar- 
ent that using Rf le > 15 MW was the better solution (com- 
pared to 20 MW). The fourth digit setting is used for con- 
trolling the number of MODIS fire pixels within the center 
aerosol pixel (F_pcounts), with “1” and “2” designating one 
and two-or-more fire pixels, respectively. It was noted that 
setting F_pcounts > 2 seems to produce similar effects on 
C e scatter plots as setting the minimum FRP value because 
both tend to eliminate small fires that potentially have under- 
estimated FRP values. The fifth digit corresponds to thresh- 
olds imposed on fire pixel counts like the fourth digit ex- 
cept that it refers to surrounding aerosol pixels in the 3x3 
aerosol-pixel matrix other than the central one. Two param- 
eters are used: setting “1” counts all the fire pixels within 
all eight aerosol pixels immediately surrounding the cen- 
tral one (F_pcounts_nearby), and setting “2” counts all the 


fires within the downwind pixels excluding the central one 
(F_pcounts_DW3). This last setting was studied as a pos- 
sible method to ensure that there is no background aerosol 
contamination from spurious plume sources that are not well 
dispersed within the 3x3 aerosol-pixel matrix. 

Table 3 shows a summary of the overall sensitivity of each 
parameter to the various threshold settings in Table 2. The 
analysis was based on global MODIS-Aqua retrievals for the 
first day of each month in 2010, for which the total num- 
ber of retrievals over this data set without any filtering was 
43 211, whereas the number of valid retrievals (after apply- 
ing the 00000 filter to ensure that valid values exist for both 
R sd and /Arc) was 28 494. Thus, roughly 34 % of the recorded 
data is invalid. The values in the table are the percentages of 
the data remaining after applying each of the thresholds. It is 
evident that the amount of available data severely decreases 
as more and more restrictions are applied. Therefore, a much 
more detailed analysis was required in order to determine the 
best choice in settings to use in the final product. After a care- 
ful evaluation of the different filters, considering their effects 
on the point scatter on plots of R sa against Rf re and the as- 
sociated correlations of the linear regression fitting vis-a-vis 
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Table 2. Value ranges of the threshold parameters in Table 1 and the combinations of their threshold settings used to derive the different 
five-digit filter configurations (00000, 10000, 1 1000, etc.) that were applied in screening out potentially erroneous or corrupted data during 
the grid-level data analysis described in Sect. 4.3. 


Parameter 

Range 

00000 

10000 

11000 

11100 

11200 

11300 

11400 

11310 

11320 

11321 

11322 

A scan angle 

[0,55] 


<30 

<30 

<30 

<30 

<30 

<30 

<30 

<30 

<30 

<30 

M wind speed 

[0,oo] 


>2 

>2 

>2 

>2 

>2 

>2 

>2 

>2 

>2 

>2 

A retrievals nearby 

[0,8] 


= 8 

= 8 

= 8 

= 8 

= 8 

= 8 

= 8 

= 8 

= 8 

= 8 

A AOT550 fire 

[0,5.05*] 




> .01 

>.01 

> .02 

>.02 

>.02 

>.02 

>.02 

>.02 

A QA AOT total 

[0,3] 


> 1 

> 1 

> 1 

> 1 

> 1 

> 1 

> 1 

> 1 

> 1 

> 1 

A_QA AOT bkgd 

[0,3] 


> 1 

> 1 

> 1 

> 1 

> 1 

> 1 

> 1 

> 1 

> 1 

> 1 

A AOT550 retr total 

[0,4] 


= 4 

= 4 

= 4 

= 4 

= 4 

= 4 

= 4 

= 4 

= 4 

= 4 

A AOT550 retr bkgd 

[0,5] 


= 5 

= 5 

= 5 

= 5 

= 5 

= 5 

= 5 

= 5 

= 5 

= 5 

A cloud frac mean 

[0,100] 



= 0 

= 0 

= 0 

= 0 

= 0 

= 0 

= 0 

= 0 

= 0 

F_pcounts 

[1,100] 








= 1 

>2 

>2 

>2 

F_pcounts_nearby 

[0,800] 










= 0 


F_pcounts DW3 

[0,300] 











= 0 

Fpower 

[0,sa 20 500*] 

>0 

>0 

>0 

> 15 

>20 

> 15 

>20 

> 15 

> 15 

> 15 

> 15 

Rsa 

[0,oo] 

>0 

>0 

>0 

>0 

>0 

>0 

>0 

>0 

>0 

>0 

>0 


* Value estimated from computations based on sensor specifications and observation geometry. 


the percentage of available valid data, 11300 was selected 
for generating the final C e product. Table 3 reports that only 
about 10 % of the available valid data is used to generate C e 
with the 11300 setting, but the confidence in the resulting C e 
values is increased by a satisfactory amount while retaining 
enough data for product development. 

This systematic data filtration process in conjunction with 
the algorithmic improvements in r.| 550 calculations described 
in Sect. 4.2 have resulted in about a 67% drop globally in 
T a 550 (which directly affects R sa and C e ), as will be seen in 
Table 6. This is a significant improvement over the Ichoku 
and Kaufman (2005) method, whose C e values were found 
to be overestimated (Sofiev et al., 2009; Kaiser et al., 2012), 
as discussed in Sect. 4.6. However, these results are still sus- 
ceptible to uncertainty and bias, as fires located in the neigh- 
boring aerosol pixels are not specifically accounted for. Al- 
though a provision was made to filter out such cases, this 
was not implemented because of high data reduction without 
a significant improvement in the result. This step will be re- 
evaluated for possible implementation in future versions of 
the FEER C e algorithm. 

4.4 Third stage: generation of smoke emission 
coefficients 

Scatter plots of R sa against Rf ie were generated for each 
1° x 1° grid cell, as illustrated in Fig. 3, using all available 
MOD1S data for the period of 2003-2010 after filtering as 
described in Sect. 4.3. Scatter plots with fewer than six data 
points were discarded. A linear least-squares regression line 
passing through the origin was fitted to each scatter plot, and 
the slope and coefficient of determination (r 2 ) were calcu- 
lated. The slope is the C e value for that grid cell. However, 
the general equation of r 2 for a regular linear least-squares 
regression analysis cannot be used for this zero-intercept fit- 


ting approach. Instead, going back to the derivation of r 2 and 
making the correct adjustments, the appropriate equation de- 
scribed in Eisenhauer (2003) was used for our situation. 

Although the process of using thresholds to remove inac- 
curate data as described in Sect. 4.3 has been successful at 
retaining the relatively higher quality R sa and R r re data se- 
ries for derivation of reliable C e , in some cases there remain 
examples where a few erroneous data points that are not suc- 
cessfully detected and filtered out can constitute outliers and 
cause large errors in this process (e.g., Fig. 3b). Such out- 
liers potentially originate from undetected errors in the data 
source, such as when the existence of clouds is undetected by 
the cloud detection algorithm. In the Fig. 3b example, when 
contrasted with Fig. 3a, only one outlier out of a total of 18 
data points cause r 2 to be as low as 0. 16, and C e to be lower 
than the expected value by a factor of six. Although the effect 
of removing outliers is usually not as drastic as this example, 
the importance of applying a filter in order to remove out- 
liers from these scatter plots before generating the final C e 
product is evident. 

The process of identifying a robust outlier removal al- 
gorithm proved to be non-trivial. Regression analysis as- 
sumes linearity, independence, homoscedasticity, and nor- 
mality. Residual plots produced from data similar to those of 
Fig. 3 show violations of at least one of these requirements, 
the most persistent being the non-normality of R s a versus 
Rf re scatter plots due to the persistent positive skewness of 
the residuals. This characteristic seems to render most if not 
all mainstream outlier algorithms unusable for the current 
study. Wisnowski (2001) describes a few highly respected 
multiple outlier detection algorithms, some of which were 
tested and found to produce many false alarms with our R sa 
vs Rf re scatter plots. It became increasingly apparent that a 
custom outlier algorithm would have to be developed specif- 
ically for these data sets. A detailed empirical study was 
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Table 3. Percentages of all available data that meet the threshold requirements in Table 2. These numbers were derived using the global 
coverage of MODIS-Aqua retrievals for the first day of each month in 2010. The number of retrievals over this data set totaled 43 211, 
whereas the number of valid retrievals (where “F _power” and “Rsa” are both greater than zero, see setting 00000 in Table 2) totaled 28 494. 
The last row (“% of Valid”) shows the overall percentages based on the 00000 setting, which gives an estimate using only valid data. Values 
are color-coded in different shades of red < 15 %, 15 % < orange < 50 %, 50 % < yellow < 70 %, green > 70 %. 


Parameter 

o 

o 

o 

o 

o 

10 000 

11 000 

11 100 

11 200 

11 300 

11 400 

11 010 

11 020 

11 310 

11 320 

11 301 

11 302 

A_scan_angle 


60.4 

60.4 

60.4 

60.4 

60.4 

60.4 

60.4 

60.4 

60.4 

60.4 

60.4 

60.4 

M_wind_speed 


84.7 

84.7 

84.7 

84.7 

84.7 

84.7 

84.7 

84.7 

84.7 

84.7 

84.7 

84.7 

A_retrievals_ 

nearby 


98.0 

98.0 

98.0 

98.0 

98.0 

98.0 

98.0 

98.0 

98.0 

98.0 

98.0 

98.0 

A_AOT550_fire 




47.6 

47.6 

34.7 

34.7 



34.7 

34.7 

34.7 

34.7 

A_QA_AOT_total 


70.3 

70.3 

70.3 

70.3 

70.3 

70.3 

70.3 

70.3 

70.3 

70.3 

70.3 

70.3 

A_QA_AOT_bkgd 


70.3 

70.3 

70.3 

70.3 

70.3 

70.3 

70.3 

70.3 

70.3 

70.3 

70.3 

70.3 

A_AOT550_retr_ 

total 


62.6 

62.6 

62.6 

62.6 

62.6 

62.6 

62.6 

62.6 

62.6 

62.6 

62.6 

62.6 

A_AOT550_retr_ 

bkgd 


57.8 

57.8 

57.8 

57.8 

57.8 

57.8 

57.8 

57.8 

57.8 

57.8 

57.8 

57.8 

A_cloud_fraction 

_mean 



52.1 

52.1 

52.1 

52.1 

52.1 

52.1 

52.1 

52.1 

52.1 

52.1 

52.1 

F_pcounts 








59.8 

40.2 

59.8 

40.2 



F_pcounts_ 

nearby 












24.9 


F_pcounts_DW3 













47.5 

F_power 

100 . 

100 . 

100 . 

81.0 

71.2 

81.0 

71.2 

100 . 

100 . 

81.0 

81.0 

81.0 

81.0 

Rsa 

65.9 

65.9 

65.9 

65.9 

65.9 

65.9 

65.9 

65.9 

65.9 

65.9 

65.9 

65.9 

65.9 

% of Total 

65.9 

23.2 

17.5 

8.8 

7.5 

6.2 

5.3 

9.3 

8.3 

1.8 

4.4 

0.8 

2.1 

% of Valid 

100 

35.2 

26.6 

13.4 

11.3 

9.5 

8.1 

14.1 

12.5 

2.8 

6.7 

1.2 

3.2 


undertaken to fully understand the variety of point distribu- 
tions that can occur in our data sets and their potential im- 
pacts on C e and r 2 resulting from the linear regression fitting 
in order to develop a robust outlier removal algorithm that 
would be optimal for our data set. The central idea behind 
the resulting outlier algorithm is to compare the fraction of 
mean squared error (MSE) measurements between the scatter 
plot with all points and without potential outliers against an 
empirically developed function in order to properly identify 
outliers. This outlier algorithm was then applied to 110 test 
scatter plots, each of which was manually assigned to one of 
15 identified scatter-point distribution categories, in order to 
rate its performance. Overall, outliers were correctly identi- 
fied and removed in 75% of the 110 cases tested, although 
three of the 15 types of scatter-point distributions showed 


a high failure rate. However, the fact that 75 % of available 
linear regression lines with outlier contamination can be rec- 
tified using this algorithm is still a vast improvement over the 
conventional outlier removal algorithms that were tested. 

When this outlier algorithm is applied to the full data set 
from both Terra and Aqua, the outlier detection rate is very 
consistent at around 30 %, regardless of the filter setting (as 
described in Sect. 4.3) that is used. If these outliers are cor- 
rectly identified, then combined with the earlier conclusion 
that 75 % of contaminated grid cells are identified by the al- 
gorithm, it is deduced that roughly 40 % of all grid cells con- 
tain outliers. Figure 6 offers an informative display of how 
the application of this outlier algorithm impacts the final C e 
product. After outlier removal, the distribution of C e values 
shifts noticeably towards higher values. This would be the 
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Figure 6. Percentages of binned C e values from all scatter plots 
based on Aqua-MODIS data and the 11300 filter (see Tables 1, 2, 
and 3) before and after outlier removal. C e appears to have shifted 
towards higher values overall. 

expected behavior for successful outlier detection since out- 
liers below the regression line (and close to the independent 
axis) have a very significant influence on the linear least- 
squares fit as compared to outliers above the line and close to 
the dependent axis. 

Initially, the scatter plots and associated linear regression 
fitting and calculations were done separately for Terra and 
Aqua data. Although a majority of the plots showed agree- 
ment between Terra and Aqua, we decided to combine them 
for deriving the final C e product. This combination offered 
two advantages: (1) it increased the number of data points 
on scatter plots with an insufficient amount of data due to 
the filtering performed above (Sect. 4.3) such that C e val- 
ues could be determined, and (2) it avoided the necessity to 
develop methods of conducting weighted averaging between 
two independent C e values for each grid cell. The resulting 
C e product is shown in Fig. 7 along with its corresponding 
r 2 map. 

4.5 Gap Ailing and quality assurance 

This polished C e product presented in Sect. 4.4 and Fig. 7 
offers the advantage of including only the highest confidence 
data, since it is based on the stringent 1 1300 filter and outlier- 
removal processes. Flowever, the tight constraints imposed 
by these processes have the effect of limiting the data suitable 
for the final product generation, such that many parts of the 
world that are known to be affected by fire do not have C e 
values generated, despite the efforts to increase coverage by 
combining Terra and Aqua data into one input stream. The 
concern of having incomplete coverage is that if a significant 
fire event were to occur in an important region, it may not 
be possible to make even a rough estimation of the smoke 
emission rates. Therefore, it is evident that some sort of filled 
product is needed. 

The possibility of a gap-filled product whereby missing C e 
values would be determined by interpolation using surround- 
ing existing values for similar land cover types was initially 


pursued. However, this procedure could not be applied at first 
because the gaps are quite extensive in certain areas, with un- 
reasonably great distances between the grid cells that need to 
be filled and those containing valid data from which their 
values can be interpolated. Thus, gaps were first filled in, 
as much as possible, using C e values based on successively 
lower filter settings starting from the 11300 setting (see Ta- 
bles 2 and 3) such that those with higher quality but less data 
are utilized before moving to those with lower quality and 
more data. To account for the differences in quality intro- 
duced by this procedure, a quality assurance (QA) product is 
provided in conjunction with the filled C e product, to serve 
as an indication of its reliability as well as to give users flex- 
ibility in the application of this product. 

The compilation process begins with the 11300-filter- 
based C e product, which is the highest confidence product, 
and progressively fills in missing data with products of lesser 
confidence: first 11000, then 10000 and finally 00000. The 
outlier removal algorithm has been applied to all except the 
00000 product. A QA flag of 0 is assigned to the lowest con- 
fidence product (i.e., 00000) and steps up to 3 for the high- 
est confidence product (i.e., 11300). C e values of the 11300 
product with r 2 > 0.7 are assigned the highest QA value of 4. 
In order to account for cases of low data availability during 
this filling process, grid cells that are already filled may be 
replaced with values from the lesser confidence product un- 
der certain conditions. The decision to replace such existing 
values is determined based on the number of data points used 
to determine C e for the previously filled value, /\'r, and for the 
new value, N n , and based on their respective r 2 values, such 
that the conditions, 

N f < AW; N n > Nr, r 2 > r 2 (14) 

must all be met, where Aii m it represents the minimum num- 
ber of data points needed to confidently fit a linear regression 
line, set to 30, which is the conventional minimum sample 
size for statistical significance. It is pertinent to recall that 
any scatter plot with less than the bare minimum of six data 
points is discarded. If Eq. (14) is satisfied for a given grid 
cell, then the C e value in the current grid cell of the new (less 
filtered) product is substituted for the existing value in the 
filled product. Likewise, the QA of the filled data is replaced 
with that of the new data. 

Finally, as many of the gaps remaining in the filled prod- 
uct as possible are filled using the C e values in nearby grid 
cells with identical land cover types. Land cover type may 
vary significantly within a grid cell at the spatial resolution 
of 1° x 1° used in this product, which can cause issues es- 
pecially since the dominant land cover type within a given 
grid cell may very well not be the one that burns most of- 
ten. Thus, the MODIS ecosystem classification map for 2004 
at 1 arcmin resolution was used to develop a custom land 
cover product at 1° x 1° resolution that reports the dominant 
fire-prone land cover type, which is used in the following 
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Figure 7. (a) The coefficient of emission (C e ) product based on MODIS 2003-2010 FRP and AOT observations from Terra and Aqua, after 
applying the 11 300 filter setting (Table 2) and outlier removal processing steps described in Sects. 4.3 and 4.4., respectively, and (b) the 
corresponding coefficient of determination (r 2 ) map. 


analysis. Grid cells that are potentially vegetated (not com- 
pletely classified as water, barren, or snow/ice) are identified 
as candidates for gap-filling, and carefully analyzed. First, a 
15x15 grid cell box is drawn around each candidate grid 
cell, and all grid cells with valid C e within this box and 
whose fire-prone land cover type is identical to that of the 
center grid cell are selected. The QA values of these selected 
grid cells are observed, and the highest minimum QA value 
(QAmin) is set such that there will be at least eight total qual- 
ified grid cells with QA > QA m i n . If this condition cannot be 
met, then no gap-filling procedure is completed in that case. 
This QA requirement is a method of balancing quantity with 
quality of data to get the most certainty in the results. Then, 
starting with a 3 x 3 box centered around the grid cell whose 
C e is being filled, the box is expanded only as necessary 
(up to a maximum of 15 x 15 window) until it contains at 
least eight grid cells that have valid C e values, the same fire- 
prone land cover type, and QA values greater than or equal to 
QA m in. If there are a sufficient number of grid cells that meet 
these criteria, the C e values of these grid cells are averaged 
and used to fill in the missing value. The gap-filled grid cell 
is assigned a QA value of zero, irrespective of those of the 
source grid cells. 

The final global 1° x 1° gridded C e product (Fig. 8a) has 
much better spatial coverage than the original (Fig. 7). The 
land areas that are not covered seem to comprise only desert 
and snow/ice regions, except for the farthest reaches of east- 
ern Russia where the last gap-filling procedure did not have a 
large enough extent to fill that area. Nevertheless, this prod- 
uct provides sufficient coverage for nearly 100 % of all veg- 
etation fires that might occur around the globe. Furthermore, 
the corresponding QA and r 2 products (Fig. 8b, c) provide 
the user with parameters for determining how accurate the 
C e in a given area might be. Table 4 shows the relative QA 
and r 2 distribution of all 13919 grid cells with C e values. 
With the exception of the grid cells whose QA equals zero, 
the majority of grid cells within each QA category are either 
in the 0.7-1 or 0.5-0. 7 r 2 range. Therefore, a user can apply 
the QA (and r 2 ) values as a filter to select only the C e val- 
ues that meet or exceed the minimum quality requirement for 


Table 4. The distribution of FEER.vl C e product QA flags among 
different coefficients of determination (r 2 ) ranges. The “N/A” cat- 
egory is for the 1° x 1° grid cells that have been filled in using the 
gap-filling process due to lack of sufficient data (and therefore lack 
an r 2 value). The r 2 range containing the most grid cells is shown 
in boldface type for each QA flag. 



r 2 


QA 

0.7-1 

0.5 -0.7 

0.3 -0.5 

0-0.3 

N/A 

Total 

4 

1745 

0 

0 

0 

- 

1745 

3 

0 

1646 

773 

182 

~ 

2601 

2 

548 

420 

321 

155 

~ 

1444 

1 

378 

395 

291 

119 

~ 

1183 

0 

243 

470 

921 

1556 

3756 

6946 

Total 

2914 

2931 

2306 

2012 

3756 

13 919 


specific applications. On the other hand, if a major fire occurs 
in a grid cell with a QA value of zero, an emissions estimate 
can still be derived, as long as the user recognizes that it is in 
fact a low-confidence estimate. This C e product is being re- 
leased as the Fire Energetics and Emissions Research version 
1 (FEER.vl) product. 

4.6 Relative evaluation of the FEER.vl C e product 

Although there is no equivalent “ground truth” data to 
validate the new FEER.vl gridded C e product, the lat- 
ter still requires a certain level of evaluation to determine 
where it stands in the spectrum of existing comparable 
products/parameters. This was done by comparing the new 
FEER.vl C e data to regional values of C e that were reported 
for 19 different regions in Table II of Ichoku and Kaufman 
(2005), hereafter referred to as “IK05”. Since the FEER C e 
product is gridded at 1° x 1°, it became necessary to generate 
a comparative set of average C e values that fit the 19 regions 
for comparison against the IK05 values. Simply averaging 
the C e grid cells within each region is unrealistic due to the 
fact that the spatial distribution of fires within each region is 
non-uniform and the certainty of the C e varies. Therefore, a 
weighted average of C e based on the number of fires within 


Atmos. Chem. Phys., 14, 6643-6667, 2014 


www.atmos-chem-phys.net/14/6643/2014/ 


C. Ichoku and L. Ellison: Global top-down smoke-aerosol emissions estimation 


6657 



o.io 

0.09 

0.08 

0.07 


0.06 : 


0.03 

U 1 0.02 
0.01 
0.00 



Figure 8. (a) The gap-filled, combined Terra and Aqua, global 1° x 1° coefficient of emission (C e ) product along with (b) the corresponding 
quality assurance (QA) map and (c) the coefficient of determination (r 2 ) map. These products were based on MODIS 2003-2010 observations 
of FRP and AOT from Terra and Aqua. 


each grid cell and also on the QA was used to generate the 
mean and standard deviation of the C e values within each 
region. 

Table 5 shows that the C e average values from the 
FEER.vl product are distinctly lower than those of 1K05 by 
a factor of 2-4.5, with the exception of East Kazakhstan, 
where they are practically equal. It is pertinent to mention 
that Ichoku and Kaufman (2005) estimated that C e values 
were probably overestimated by a factor of 2, and Sofiev 
et al. (2009) by applying a more rigorous plume dispersion 
modeling found C e values that were lower than those of IK05 
by a factor of 2 to 3. Kaiser et al (2012) also found values 
that were lower than those of IK05. The fact that those sub- 
sequent studies, including the current study, produced lower 
values than those of 1K05, confirms that 1K05 values were 
indeed probably overestimated and suggests that those from 
the current study are more realistic. The change from 1K05 to 
the current study can be categorized into two types, namely, 
algorithms and input data versions/sources. It is necessary 
to characterize these two types of change independently in 
order to determine their relative contributions (as will be re- 
ported in Table 6). 

To account for the effects of using new input data ver- 
sions/sources, an updated version of the IK05 product (here- 
after referred to as IKu) was generated by ingesting the 
new data being used in the current study into an algorithm 
that matches that of 1K05 as closely as possible. Recall that 
the 1K05 C e values were based on the MODIS collection 4 
FRP and AOT products, with wind data from the NCEP re- 
analysis data set (GDAS1). By contrast, the IKu C e values 
are based on the MODIS collection 5 FRP and AOT prod- 
ucts, with wind data from the MERRA reanalysis data set. 
Differences in C e from IK05 to IKu should only be due to 
changes in data versions and sources, whereas the effects of 
the algorithmic alterations described in Sect. 4.2 can be iso- 
lated by comparing IKu to FEER.vl. 


Using the relationships defined in Sect. 4.1, it is evident 
that 


C e oc 


Rsa 

FRP 

Afca 


oc 


T FRP 


oc 


M sa • WS 
L FRP 


oc 


Mi- A- WS 
L FRP 


r f - A-WS 
L ■ FRP 


(15) 


In other words, C e is directly proportional to the fire-emitted 
AOT (r,[ 55() ), aerosol-pixel area and wind speed, but inversely 
proportional to the plume length and FRP. Three of the 
five variables on the right-hand side of Eq. (15) (r a r 55() , WS 
and FRP) have updated data sources in FEER.vl, and three 
(r.[ 550 , A and L ) have updated derivations. However, both 
A and L, which are dependent on each other, can be ad- 
justed together here to emulate the IK05 algorithm such that 
A would be equal to the area of only one aerosol pixel, and 
L would be halved (using IK05 definitions, Tikos / Tfeer = 


/4 A — 0.5). These adjustments are made for the fol- 


lowing analysis. If the ratio of a variable in FEER.vl to the 
same in IK05 is represented by R, then 


*Ce = 


Ce,FEER 

C e ,IK05 


( Tf-AWS \ 

V LFRP /feer 
/ r f -AWS \ 

V ™ Akos 


Rx f ■ Rws 

Rl ■ Rfrp 


(16) 


Equation (16) quantifies the change in C e due to both input 
source and algorithmic alterations from IK05 to FEER.vl. 
Changes in only data sources from IK05 to IKu are captured 
in the relationship: 


Rc e = 


Rtf • Rws 

Rfrp 


(17) 


because the calculation of L does not involve the use of data 
from different sources. The relationship that quantifies the 
algorithmic changes from IKu to FEER.vl is given by 


R 


Tf 


R Ce = — 
e Rl 


(18) 


www.atmos-chem-phys.net/14/6643/2014/ 


Atmos. Chem. Phys., 14, 6643-6667, 2014 




6658 


C. Ichoku and L. Ellison: Global top-down smoke-aerosol emissions estimation 


Table 5. Estimates of regional FRE-based smoke-aerosol emission coefficients (C e ) from MODIS are shown here for different regions using 
both the original method as reported in Ichoku and Kaufman (2005) and version 1 of the new FEER C e product (FEER.vl). 




PM emission coefficients 

Region 

Description 


(kg M J _ 1 ) 




IK05 

FEER.vl 



calculated mean 

st.dev. 

Savanna and grassland regions 

Brazil-Cer 

Brazil Cerrado savanna region 

0.048 

0.016 

0.009 

S. America 

South America below 20° S 

0.061 

0.020 

0.013 

W. Africa 

West Africa 

0.059 

0.021 

0.012 

Zambia 

Zambia in southern Africa 

0.076 

0.018 

0.005 

Tropical forest regions 

Borneo 

Borneo Island of Indonesia 

0.079 

0.032 

0.019 

Brazil-For 

Brazil tropical forest region 

0.063 

0.019 

0.009 

Celebes-Moluccas 

Celebes and Moluccas Islands, Indonesia 

0.068 

0.028 

0.020 

Congo 

Congo tropical forest, Africa 

0.048 

0.015 

0.006 

Boreal forest regions 

Alaska 

Alaska 

0.020 

0.012 

0.016 

Canada 

Canada below 70° N (excluding Quebec) 

0.020 

0.012 

0.013 

Quebec 

Quebec and eastern Ontario 

0.020 

0.009 

0.021 

Siberia 

Siberia North of 60° N 

0.057 

0.024 

0.018 

Cropland/natural vegetation regions 

Moscow 

Moscow and environs 

0.100 

0.026 

0.011 

S. Russia 

Southern Russia 

0.084 

0.018 

0.007 

St. Petersburg 

St. Petersburg and environs 

0.104 

0.023 

0.009 

Unclassified 

Europe 

Europe (excluding Russia) 

0.056 

0.024 

0.017 

E. Kazakhstan 

East Kazakhstan 

0.018 

0.019 

0.011 

Mongolia 

Mongolia 

0.033 

0.022 

0.014 

Philippines 

The Philippines 

0.127 

0.039 

0.024 


because the way in which FRP and wind speed are calculated 
remains the same between the two algorithms. 

The relationships shown in Eqs. (16), (17) and (18) can be 
utilized to test whether the differences between the IK05, IKu 
and FEER.vl product data sets can fully explain the change 
in C e between IK05 and FEER.vl shown in Table 5 as well 
as to identify the main factor responsible for the change - 
change in algorithm or the input data version/source. The 
only available data from the original IK05 data set cover rel- 
atively short time periods (Terra: 25 June 2002 to 4 October 
2002, and Aqua: 25 June 2002 to 31 December 2002). The 
fact that these ranges do not cover a full year means that any 
seasonal differences that may exist will be lost and will there- 
fore cause the resulting data to be biased low or high. Nev- 
ertheless, these 2002 data sets were used to estimate K( c , 
by first pairing corresponding individual data points in the 
IK05, IKu and FEER.vl data sets. The ratios between IK05 
and IKu, between IKu and FEER.vl, and between IK05 and 


FEER.vl were calculated for each data point for AOT, wind 
speed, FRP, and plume length. Subsequently, the ratio of C e 
was calculated for each data point pair according to Eq. (16) 
for the transition from IK05 to FEER.vl, Eq. (17) for the 
transition from IK05 to IKu, and Eq. (18) for the transi- 
tion from IKu to FEER.vl. To appropriately represent these 
matched data points and ratios in a uniform fashion within 
the spatial domains outlined in Table 5, they were binned 
into a global grid at a spatial resolution of 0.5 x 0.5° and 
then filtered according to the appropriate settings reported in 
Table 2 using the QA values from the FEER.vl C e product 
in Fig. 8. Finally, the median of those ratios within each grid 
cell was retrieved, and the mean of these values (weighted 
identically as was done in Table 5 for C e ) within each region 
were reported, as displayed in Table 6. 

In Table 6, column 1 (highlighted yellow) shows ob- 
served changes in C e from IK05 to FEER.vl. The subse- 
quent columns outline the process of deriving the predicted 
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Table 6. The observed changes of C e and intermediate parameters from the original Ichoku and Kaufman (2005) method (IK05) to FEER.vl 
is shown here for the regions listed in Table 5. The values in the yellow highlighted column on the left-hand side are the observed changes 
in C e from IK05 to FEER.vl from Table 5. The subsequent columns outline the sample size and mean parameter changes during the process 
of deriving the predicted changes in C e from IK05 to FEER.vl according to Eqs. (16), (17), and (18), the results of which are highlighted in 
orange and yellow. The “IK=> FEER.vl” on the first header row indicates that both the IKu to FEER.vl and IK05 to FEER.vl transitions 
are included in the columns beneath. Both Terra and Aqua data were used in these calculations. 




IK05 => IKu 


IK => FEER.vl 


Region 

Ce 

(FEER/ 

IK05) 

N FRP WS AOT 

Ce 

(IKu/ 

IK05) 

AOT AOT 
(FEER (FEER 
N FRP WS L /IKu) /IK05) 

Ce Ce 

(FEER (FEER 
/IKu) /IK05) 


Savanna and grassland regions 


Brazil-Cer 

0.33 

851 

1.61 

0.88 

0.99 

0.56 

295 

1.48 

0.97 

0.85 

0.35 

0.38 

0.41 

0.26 

S. America 

0.32 

897 

1.57 

0.86 

0.99 

0.55 

385 

1.44 

0.91 

0.79 

0.32 

0.31 

0.41 

0.23 

W. Africa 

0.35 

766 

1.60 

1.02 

0.96 

0.62 

565 

1.50 

1.03 

0.82 

0.35 

0.32 

0.43 

0.25 

Zambia 

0.23 

518 

1.59 

0.79 

0.99 

0.50 

470 

1.46 

0.90 

0.85 

0.35 

0.32 

0.41 

0.19 


Tropical forest regions 


Borneo 

0.40 

193 

1.34 

1.03 

1.02 

0.82 

98 

1.26 

1.03 

0.83 

0.29 

0.35 

0.35 

0.37 

Brazil-For 

0.29 

832 

1.50 

0.72 

1.04 

0.51 

391 

1.50 

0.83 

0.83 

0.31 

0.36 

0.38 

0.21 

Celebes- 

Moluccas 

0.41 

116 

1.28 

0.66 

0.99 

0.52 

81 

1.25 

0.86 

0.81 

0.30 

0.30 

0.38 

0.25 

Congo 

0.31 

913 

1.65 

0.91 

1.00 

0.56 

621 

1.50 

1.04 

0.84 

0.33 

0.32 

0.40 

0.21 

Boreal forest 

regions 














Alaska 

0.62 

26 

2.32 

1.14 

1.00 

0.55 

12 

2.01 

1.14 

0.81 

0.29 

0.25 

0.36 

0.21 

Canada 

0.58 

65 

2.39 

0.99 

0.97 

0.46 

31 

2.36 

0.93 

0.82 

0.33 

0.35 

0.41 

0.20 

Quebec 

0.46 

54 

2.10 

0.82 

0.98 

0.43 

35 

1.83 

0.85 

0.82 

0.28 

0.23 

0.34 

0.14 

Siberia 

0.42 

290 

1.65 

0.86 

0.98 

0.53 

174 

1.63 

0.87 

0.82 

0.32 

0.27 

0.39 

0.18 


Cropland/natural vegetation regions 


Moscow 

0.26 

123 

1.38 

1.01 

1.00 

0.76 

60 

1.32 

1.09 

0.85 

0.29 

0.29 

0.34 

S. Russia 

0.21 

101 

1.59 

0.94 

0.97 

0.60 

62 

1.42 

1.03 

0.85 

0.34 

0.29 

0.40 

St. Petersburg 

0.23 

58 

1.33 

1.02 

1.01 

0.78 

35 

1.27 

1.20 

0.84 

0.30 

0.30 

0.35 


Unclassified 


Europe 
E. Kazakhstan 
Mongolia 
Philippines 

0.43 

1.06 

0.67 

0.30 

544 1.49 0.94 0.96 
431 1.72 1.01 0.86 
149 1.59 0.91 0.88 
20 1.17 0.86 0.89 

0.64 

0.55 

0.54 

0.68 

172 1.40 0.97 0.84 0.29 0.29 

152 1.50 1.09 0.84 0.30 0.27 

46 1.62 0.96 0.85 0.31 0.31 

6 1.06 0.87 0.79 0.23 0.19 

0.35 

0.37 

0.36 

0.29 

0.24 

0.21 

0.24 

0.22 







Global avg. 


13312 1.62 0.90 0.98 

0.56 

6452 1.48 1.00 0.83 0.33 0.32 

0.40 0.22 


changes in C e from 1K05 to FEER.vl according to Eqs. 
(16), (17), and (18), the results of which are shown in 
the last column (highlighted yellow). Both Terra and Aqua 
data were used in these calculations. The two main pro- 
cess changes have been separated out: columns 2-6 (labeled 
“IK05 =>■ IKu”) clearly showing the effect of altering only the 
data inputs, and columns 7-14 (labeled “IK =£■ FEER.vl”) 
showing both the effect of altering only the algorithm 
(IKu =>• FEER.vl) and the combination of algorithm alter- 
ation and data updates (IK05 =>■ FEER.vl). From the result- 
ing maps showing the global variation in ratios of r,|- 5(| , WS, 
FRP, and L. it was apparent that the change in each variable 
is uniform throughout the globe. On average, the change in 
C e due to differing data sources is about a 40% decrease 
(column 6), mostly from the change in FRP from collec- 
tion 4 to 5 (i.e., without and with multiplication by fire -pixel 
area, respectively), whereas the algorithm alterations cause 
about a 60% decrease in C e (column 13), resulting in an 


overall decrease in C e of about 80% globally (column 14). 
Even though these combined effects of data-source and al- 
gorithm changes are slightly overcompensating compared to 
the observed differences listed in column 1, it can be stated 
that the observed reduction in C e values between the IK05 
and FEER.vl is indeed realistic. The changes in wind speed, 
plume distance and FRP due to algorithmic changes are small 
relative to the large change in AOT. Therefore, most of the 
change in C e is attributable to the change in fire-emitted 
AOT (r a f 550 ). Figure 9 shows the global distribution of r a550 
changes due to data version/source change (i.e., from IK05 
to IKu, Fig. 9a) and due to algorithm change (i.e., from IKu 
to FEER.vl, Fig. 9b). Interestingly, when the new collection 
5 AOT data are used in lieu of collection 4, r a550 actually in- 
creases in most cases around the globe, confirming that the 
lower C e values from IK05 to FEER.vl due to AOT is very 
strongly attributable to the change in the r.[ 550 algorithm. In 
fact, the ratios of t a550 in the data-source part are very near 
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Figure 9. The ratio of fire-emitted AOT values at 550 nm wavelength (t^q) between the new (FEER.vl) and old (IK05) products mapped 
on a 0.5 x 0.5° global grid. The changes in are due only to (a) upgrading the data source from collection 4 to collection 5, and 
(b) algorithmic changes. 


unity (column 5, Fig. 9a), whereas in the algorithm alteration 
part it is around 0.3 (columns 11 and 12, Fig. 9b). It is per- 
tinent to recall that the r,| 550 algorithmic changes mainly in- 
volve (1) using wind direction to detennine which AOT val- 
ues to classify as plume or background, and (2) taking the 
average of the upwind AOT values (instead of just the min- 
imum value) as the background in an effort to account for 
contamination from external aerosols. Although these modi- 
fications have resulted in a severe change in the derived r.| 55() , 
this change translates to increased confidence in C e . 

5 Emissions calculation results 

The new FEER.vl C e product is used to demonstrate the top- 
down derivation of emission rates and totals from satellite 
measurements of FRP. The resulting emissions are compared 
against other emission inventories to gain a general under- 
standing for how model simulations will change when using 
this new FEER.vl inventory. 

5.1 Emissions estimates (rates and totals) 

The FEER.vl C e product is used to derive smoke-aerosol 
emissions by simple multiplication, as represented in Eq. (2) 
and the associated discussion. When C e (kgMJ -1 ) is mul- 
tiplied directly by FRP (in MW or MJs -1 ), instantaneous 
emission rates (in kgs -1 ) are derived, whereas when mul- 
tiplied by FRE (in MJ) representing a finite (e.g., daily, 
monthly, or yearly) time period, the result is emission totals 
(in kg) for that time period. Generating a global FRE prod- 
uct for use in this analysis is not straightforward due to the 
fact that semi-continuous measurements of unsaturated FRP 
around the entire globe is not currently available, though it 
is expected that this situation will improve within the next 
decade or so, given the anticipated launches of different geo- 
stationary and polar-orbiting satellite missions by some of 
the major space agencies. However, to closely compare emis- 
sions based on the new FEER.vl C e product with other emis- 
sions products, this study uses FRP data from the 0.5° x 0.5° 
gridded monthly data set derived from MODIS observations 


aboard the Terra and Aqua satellites as part of the GFAS.vl 
product (http://gmes-atmosphere.eu/fire, Kaiser et al., 2012). 

The GFAS.vl values of monthly average FRP in W m -2 
were simply multiplied by the number of days in each calen- 
dar month to get FRE in J m -2 , as was done in the GFAS.vl 
algorithm (Kaiser et ah, 2012). Such derivation of monthly 
average FRE based on only four or less MODIS fire obser- 
vations a day (from Terra and Aqua satellites) result in high 
uncertainty, as it cannot capture the fire diurnal cycle. How- 
ever, that is currently the only feasible way to obtain FRE 
globally. Higher- frequency (sub-hourly) observations from 
a few available geostationary satellite sensors that measure 
FRP have different characteristics and produce an average of 
17-38 % underestimation relative to MODIS coincident FRP 
observations (Roberts et ah, 2005; Xu et ah, 2010). More- 
over, a combination of these geostationary FRP data still 
does not provide global coverage, as some large biomass- 
burning regions, including Siberia, Central Asia, and India, 
are left uncovered (Zhang et ah, 2012). Since the GFAS.vl- 
based FRE data are global, publicly available, and being used 
in the European Union’s Monitoring Atmospheric Composi- 
tion and Climate (MACC) project (http://gmes-atmosphere. 
eu/fire), they were considered appropriate for use in de- 
riving emissions using the FEER.vl C e product to enable 
comparison with existing emissions inventories, as described 
below. Therefore, these GFAS.vl monthly FRE values at 
0.5° x 0.5° resolution were multiplied by the FEER.vl C e 
product at 1° x 1° resolution to obtain the monthly emissions 
of smoke aerosols around the globe at 0.5° x 0.5° resolution. 
Then, the monthly emissions for all months of a calendar 
year were summed up to get yearly emissions estimates, such 
as the 2010 example shown in Fig. 10. 

5.2 Comparison with other emissions inventories 

The FEER.vl monthly emissions were compared with some 
of the existing emissions products - GFED.v3, GFAS.vl, and 
QFED.v2 - as a way of evaluating the FEER.vl emissions 
within the context of these existing global emission invento- 
ries that are currently used by the research and operational 
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Total emissions: 72.8 Tg 



Mass of emission 1 g/m* ] 

Figure 10. FEER.vl emissions estimates of total particulate mat- 
ter (TPM) for all of 2010 on a 0.5° x 0.5° resolution global grid. 
These values are generated from Eq. (2), using the FEER.vl C e 
product in conjunction with FRE data derived by multiplying the 
GFAS.vl monthly average FRP data by the total number of days in 
each month. 


communities. It should be noted that there are a few dis- 
similarities between these products. First, like FEER.vl, 
QFED.v2 is based on the top-down approach using satel- 
lite measurements of both aerosols and FRP, whereas both 
GFED.v3 and GFAS.vl are based on the bottom-up approach 
using literature-extracted emission factors (EF) to multiply 
burned biomass estimates from satellite observations of FRP 
(GFAS.vl) or fire pixel counts and burned areas (GFED.v3). 
Secondly, the emissions values used for comparison from 
both GFED.v3 and GFAS.vl represent the smoke TPM emis- 
sions, whereas for QFED.v2, whose product exists as the 
component species of smoke aerosols, the closest equivalent 
product is particulate matter < 2.5 pm aerodynamic diameter 
(PM 2 . 5 ). The ratio of PM 2.5 to TPM (by ratioing their cor- 
responding emission factors) is estimated to range between 
65% and ~ 100% depending on ecosystem type (e.g., An- 
dreae and Merlet, 2001; Akagi et al., 2011). Thus, ideally, 
smoke PM 2.5 emissions (from QFED.v2) for a given area 
and time period should be expected to be lower than the cor- 
responding TPM emissions (from FEER.vl, GFED.v3, and 
GFAS.vl). These different data sets were aggregated region- 
ally according to the regional biomass burning partitions pro- 
vided in Kaiser et al. (2012) as delineated in Fig. 11, and 
plotted as time series annual total smoke TPM emissions 

(Fig. 12). 

All the emissions products portray similar temporal pat- 
terns, with lows and highs occurring in the same years, for 
both the global and regional plots (Fig. 12). This may be 
due at least in part to the fact that all products are influ- 
enced by MOD1S fire-pixel counts, either directly or indi- 
rectly. GFAS.vl emissions are generally equal to those of 
GFED.v3, because the former was scaled to match the latter 
(Kaiser et ah, 2012). Globally, GFED.v3 and GFAS.vl each 



Figure 11. Regional partitions as defined in Kaiser et al. (2012) 
that are used in this paper to compare FEER.vl emissions with 
GFED.v3, GFAS.vl, and QFED.v2 emission inventories. The back- 
ground MODIS true-color image shows fire locations (red dots) de- 
tected by MODIS from both Terra and Aqua for all of 2010, to il- 
lustrate the global spatial distribution of annual fire occurrence. 


constitutes only about 55% of the FEER.vl annual TPM 
emissions. Since the GFAS.vl FRE data set was also used for 
FEER.vl, it follows that the large difference between their 
emission products stem from the relative magnitudes of the 
C e used to generate them. Furthermore, given that it was al- 
ready established that the TPM emissions in GFAS.vl (and 
by inference also in GFED.v3) need to be boosted by a fac- 
tor of 2-4 to match realistic global distributions of aerosols, it 
follows that FEER.vl C e results are probably closer to realis- 
tic values. However, although the QFED.v2 emissions repre- 
sent only PM 2 . 5 , which should be lower than TPM, paradox- 
ically, it is slightly higher than FEER.vl global TPM emis- 
sions. 

The relationship between the FEER.vl emissions and 
those of GFED.v3, GFAS.vl, and QFED.v2 portrays signif- 
icant regional differences, as indicated by the regional plots 
in Fig. 12. In North America (NAme), incidentally, FEER.vl 
emissions seem to agree closely with those of GFED.v3 
and GFAS.vl, whereas QFED.v2 (though only PM 2 . 5 ) shows 
double the values of the former three TPM emissions in- 
ventories. Not surprisingly, out of all the regions, NAme 
has the largest distribution of the lowest QA and r 2 values 
for the FEER.vl C e values, as shown in Fig. 8 . We suspect 
that FEER.vl C e values are severely underestimated in this 
NAme region probably because, among other possible rea- 
sons, the gap-filled areas are quite extensive with very low 
QA values that may have tended toward underestimation. On 
the other hand, QFED.v2 appears to have been overestimated 
in northern and southern Asia (NAsi and SAsi), perhaps due 
to contamination from the persistent regional pollution, since 
QFED.v2 is based on regional aerosol observations in con- 
trast to FEER.vl, which is based on near-source AOT mea- 
surements. Similarly, GFED.v3 is probably overestimated in 
tropical Asia (TAsi) only in 2002 and 2006, although the 
investigation of possible reasons for these two anomalous 
years is beyond the scope of this paper. 
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Figure 12. Time series of yearly total global and regional emissions of total particulate matter (TPM) in Tg from 2000-2012 for FEER.vl, 
GFED.v3, QFED.v2 and GFAS.vl, depending on data availability for this study. QFED.v2 values (dotted line) are for PM 2 . 5 . The regions 
represented are delimited on the map in Fig. 11. 


6 Summary and conclusions 

This study has presented a first attempt at providing C e - 
an index that is similar to EF - for every 1° x 1° cell con- 
taining burnable vegetation globally. Whereas EF is used to 
multiply burned biomass estimates to calculate emissions, C e 
is the equivalent parameter used to multiply time-integrated 
satellite measurements of FRP to estimate emissions. Thus 
the FEER.vl global gridded C e product developed in this 
study for TPM emissions estimation has several important 
attributes, of which the most significant are that it (1) is 
the first global gridded product in the family of “emission 
factors”, whereas existing products specify one value per 
ecosystem type; (2) requires only direct satellite measure- 
ments of FRP or its time-integrated FRE to generate emis- 
sion rates or totals, respectively, whereas regular EF values 
still require estimation of burned biomass through an intri- 
cate process fraught with high uncertainty; and (3) is the 
only variable in the family of “emission factors” that does 
not require pre-determination of the ecosystem type of an 
actively burning fire to evaluate its emission rate in near-real 
time, which is essential for operational activities, such as the 
monitoring and forecasting of smoke emission impacts on air 
quality. 


Although the FEER.vl global gridded C e product was 
based on the original approach proposed by Ichoku and 
Kaufman (2005), this study implemented significant im- 
provements in all stages of the product development. The lat- 
est available versions (collection 5) of both the aerosol and 
fire products from MOD1S were used, along with MERRA 
meteorological data from the GEOS-5 global assimilation 
model. The identification of near-source plume and back- 
ground pixels from the MOD1S AOT data set was based 
on actual wind directions from MERRA. Rigorous methods 
were used to determine the valid ranges of all parameters uti- 
lized in the algorithm, in order to limit the effects of spuri- 
ous errors and uncertainties from measurements and assump- 
tions. These updates in data versions and algorithm resulted 
in an overall decrease in regional average C e values by a fac- 
tor of 2-4.5 relative to those of Ichoku and Kaufman (2005). 
This decrease seems reasonable, as observed by recent stud- 
ies that evaluated those C e values based on model analyses 
(e.g., Sofiev et al., 2009; Kaiser et al., 2012). Nevertheless, 
the FEER.vl global gridded C e product may still contain sev- 
eral limitations and uncertainties, which may be due to vari- 
ous factors, such as (1) uncertainties in the satellite retrievals 
of AOT and FRP, (2) omission of smaller fires or even larger 
fires that are mostly smoldering with significant smoke emis- 
sion but limited radiant energy signal below the MODIS de- 
tection limit, (3) possibility of erroneously including external 
aerosols to specific plumes being analyzed or having large 
variations in the aerosol background surrounding the plume, 
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(4) smoke underestimation due to the erroneous masking 
out of near-source thick smoke plumes as cloud during the 
aerosol retrieval process, (5) lack of knowledge of plume in- 
jection heights, ( 6 ) use of wind vectors at 850 mb globally, 
(7) uncertainty in the MERRA wind vectors used in the cal- 
culations of smoke emission rate and trajectory, ( 8 ) assump- 
tion of a single value of smoke-aerosol mass extinction ef- 
ficiency globally, and (9) uncertainties due to the gap-filling 
process of the FEER.vl global gridded C e product. There- 
fore, there is need to find ways of validating this product. 
Fortunately, the fact that this global C e product is anchored 
on a geographically fixed grid system makes validation much 
more feasible than is the case for existing EF values whose 
geographical attributes may have been lost, thereby making 
them difficult to replicate or to trace to a specific geographic 
domain. Thus, for the FEER.vl global gridded C e product, 
deliberate effort could be made to conduct field experiments 
within any 1° x 1° grid cell for use in validating its C e value. 

Pending the validation of this FEER.vl global gridded 
C e product in a representative sample of locations, perhaps 
through the use of observations in conjunction with regional 
modeling, QA flags (ranging from 0 to 4 in increasing order 
of quality) have been provided with the product to guide the 
user in using this product for different applications. These 
QA flags were based on several qualitative and quantitative 
considerations including the r 2 from the zero-intercept lin- 
ear least-squares regression fitting of smoke-aerosol emis- 
sion rates against FRP. A corresponding gridded map of r 2 
is also provided for reference. Thus, a user desiring to de- 
rive only high-quality emissions can use the QA as a filter 
to select only the C e values with the highest quality required, 
while the corresponding r 2 value can give a general idea as to 
whether this QA is based on quantitative or qualitative con- 
siderations. On the other hand, if a fire occurs in a grid cell 
for which emissions estimates are needed to determine the 
general smoke trajectory without the need for precise quan- 
titative estimates of concentrations, even C e values having a 
QA value of zero can be used to accomplish the desired task. 

The FEER.vl global gridded C e product was used to gen- 
erate global and regional total emissions of TPM, which 
were compared against existing emissions inventories, in- 
cluding GFED.v3, GFAS.vl, and QFED.v2. The FEER.vl 
annual total TPM emissions are low in North America, where 
they had comparable magnitudes as those of GFED.v3 and 
GFAS.vl, each of which was about half of the PM 2.5 emis- 
sions from QFED.v2. Pending validation, with the exception 
of the North America, FEER.vl and QFED.v2 seem com- 
parable in most regions relative to GFED.v3 and GFAS.vl 
emissions, which are considered low. Since GFED.v3 and 
GFAS.vl products are based on bottom up approaches (with 
regards to the determination of the emission factors used), 
whereas FEER.vl and QFED.v2 are based on top-down ap- 
proaches (in relation to the emission coefficients used), it 
is reasonable to assume that top-down approaches based 
on satellite measurements would yield smoke distributions 


that have a closer resemblance to satellite observations of 
aerosols. Therefore, it is recommended that increased ef- 
fort be made toward further enhancement of top-down ap- 
proaches, not only for aerosol emissions, but also for gaseous 
emissions. It is hoped that this approach will become more 
and more accurate and beneficial with continued improve- 
ment in the satellite retrievals of these aerosols and gases. 

The FEER.vl C e product, which is currently based on 
2003-2010 MOD1S observations, will certainly require fu- 
ture updates as new, improved, and more representative data 
inputs become available from other relevant sources. Further- 
more, even in places where the current C e values are reason- 
ably accurate, over time, changing land-cover conditions and 
fire regimes will invariably necessitate the revision of these 
C e values, which may indicate diurnal, seasonal, annual or 
longer-term temporal variability. Future studies will reveal 
the approaches required to ensure optimal accuracy in time 
and space. 

The current study has been focused on the development 
of a global gridded C e product for smoke TPM because it is 
based on the total columnar AOT parameter as retrieved from 
satellite observations. Although it is recognized that mod- 
eling activities often require smoke-aerosol speciation into 
its various components such as organic carbon (OC), BC, 
or PM 2 . 5 , it was beyond the scope of this study to derive 
emission coefficients for these smoke constituent species, 
as it would have involved several assumptions (with asso- 
ciated compounding uncertainties) to estimate any one of 
them from satellite AOT retrievals. However, the user of the 
FEER.vl TPM C e product may optionally estimate corre- 
sponding C e values for any of the other smoke-aerosol con- 
stituents by multiplying with their emission ratios relative to 
TPM. Such emission ratios can be obtained from the litera- 
ture or derived from the constituent emission factors, which 
are also available in the literature depending on ecosystem 
type (e.g., Andreae and Merlet, 2001; Akagi et al., 2011). 
Indeed, the FEER.vl global gridded TPM C e product devel- 
oped in this paper represents a versatile foundational product 
that can lead to several important advances in fire emissions 
research and applications. 
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